home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / integration / qc25s.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-05-05  |  5.6 KB  |  239 lines
  1. /* integration/qc25s.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. struct fn_qaws_params
  21. {
  22.   gsl_function *function;
  23.   double a;
  24.   double b;
  25.   gsl_integration_qaws_table *table;
  26. };
  27.  
  28. static double fn_qaws (double t, void *params);
  29. static double fn_qaws_L (double x, void *params);
  30. static double fn_qaws_R (double x, void *params);
  31.  
  32. static void
  33. compute_result (const double * r, const double * cheb12, const double * cheb24,
  34.         double * result12, double * result24);
  35.  
  36.  
  37. static void
  38. qc25s (gsl_function * f, double a, double b, double a1, double b1,
  39.        gsl_integration_qaws_table * t,
  40.        double *result, double *abserr, int *err_reliable);
  41.  
  42. static void
  43. qc25s (gsl_function * f, double a, double b, double a1, double b1,
  44.        gsl_integration_qaws_table * t,
  45.        double *result, double *abserr, int *err_reliable)
  46. {
  47.   gsl_function weighted_function;
  48.   struct fn_qaws_params fn_params;
  49.   
  50.   fn_params.function = f;
  51.   fn_params.a = a;
  52.   fn_params.b = b;
  53.   fn_params.table = t;
  54.  
  55.   weighted_function.params = &fn_params;
  56.     
  57.   if (a1 == a && (t->alpha != 0.0 || t->mu != 0))
  58.     {
  59.       double cheb12[13], cheb24[25];
  60.  
  61.       double factor = pow(0.5 * (b1 - a1), t->alpha + 1.0);
  62.  
  63.       weighted_function.function = &fn_qaws_R;
  64.  
  65.       gsl_integration_qcheb (&weighted_function, a1, b1, cheb12, cheb24);
  66.  
  67.       if (t->mu == 0)
  68.     {
  69.       double res12 = 0, res24 = 0;
  70.       double u = factor;
  71.  
  72.       compute_result (t->ri, cheb12, cheb24, &res12, &res24);
  73.  
  74.       *result = u * res24;
  75.       *abserr = fabs(u * (res24 - res12));
  76.     }
  77.       else 
  78.     {
  79.       double res12a = 0, res24a = 0;
  80.       double res12b = 0, res24b = 0;
  81.  
  82.       double u = factor * log(b1 - a1);
  83.       double v = factor;
  84.  
  85.       compute_result (t->ri, cheb12, cheb24, &res12a, &res24a);
  86.       compute_result (t->rg, cheb12, cheb24, &res12b, &res24b);
  87.  
  88.       *result = u * res24a + v * res24b;
  89.       *abserr = fabs(u * (res24a - res12a)) + fabs(v * (res24b - res12b));
  90.     }
  91.  
  92.       *err_reliable = 0;
  93.  
  94.       return;
  95.     }
  96.   else if (b1 == b && (t->beta != 0.0 || t->nu != 0))
  97.     {
  98.       double cheb12[13], cheb24[25];
  99.       double factor = pow(0.5 * (b1 - a1), t->beta + 1.0);
  100.  
  101.       weighted_function.function = &fn_qaws_L;
  102.  
  103.       gsl_integration_qcheb (&weighted_function, a1, b1, cheb12, cheb24);
  104.  
  105.       if (t->nu == 0)
  106.     {
  107.       double res12 = 0, res24 = 0;
  108.       double u = factor;
  109.  
  110.       compute_result (t->rj, cheb12, cheb24, &res12, &res24);
  111.  
  112.       *result = u * res24;
  113.       *abserr = fabs(u * (res24 - res12));
  114.     }
  115.       else 
  116.     {
  117.       double res12a = 0, res24a = 0;
  118.       double res12b = 0, res24b = 0;
  119.  
  120.       double u = factor * log(b1 - a1);
  121.       double v = factor;
  122.  
  123.       compute_result (t->rj, cheb12, cheb24, &res12a, &res24a);
  124.       compute_result (t->rh, cheb12, cheb24, &res12b, &res24b);
  125.  
  126.       *result = u * res24a + v * res24b;
  127.       *abserr = fabs(u * (res24a - res12a)) + fabs(v * (res24b - res12b));
  128.     }
  129.  
  130.       *err_reliable = 0;
  131.  
  132.       return;
  133.     }
  134.   else
  135.     {
  136.       double resabs, resasc;
  137.  
  138.       weighted_function.function = &fn_qaws;
  139.   
  140.       gsl_integration_qk15 (&weighted_function, a1, b1, result, abserr,
  141.                 &resabs, &resasc);
  142.  
  143.       if (*abserr == resasc)
  144.     {
  145.       *err_reliable = 0;
  146.     }
  147.       else 
  148.     {
  149.       *err_reliable = 1;
  150.     }
  151.  
  152.       return;
  153.     }
  154.  
  155. }
  156.  
  157. static double
  158. fn_qaws (double x, void *params)
  159. {
  160.   struct fn_qaws_params *p = (struct fn_qaws_params *) params;
  161.   gsl_function *f = p->function;
  162.   gsl_integration_qaws_table *t = p->table;
  163.  
  164.   double factor = 1.0;
  165.   
  166.   if (t->alpha != 0.0)
  167.     factor *= pow(x - p->a, t->alpha);
  168.  
  169.   if (t->beta != 0.0)
  170.     factor *= pow(p->b - x, t->beta);
  171.  
  172.   if (t->mu == 1)
  173.     factor *= log(x - p->a);
  174.  
  175.   if (t->nu == 1)
  176.     factor *= log(p->b - x);
  177.  
  178.   return factor * GSL_FN_EVAL (f, x);
  179. }
  180.  
  181. static double
  182. fn_qaws_L (double x, void *params)
  183. {
  184.   struct fn_qaws_params *p = (struct fn_qaws_params *) params;
  185.   gsl_function *f = p->function;
  186.   gsl_integration_qaws_table *t = p->table;
  187.  
  188.   double factor = 1.0;
  189.   
  190.   if (t->alpha != 0.0)
  191.     factor *= pow(x - p->a, t->alpha);
  192.  
  193.   if (t->mu == 1)
  194.     factor *= log(x - p->a);
  195.  
  196.   return factor * GSL_FN_EVAL (f, x);
  197. }
  198.  
  199. static double
  200. fn_qaws_R (double x, void *params)
  201. {
  202.   struct fn_qaws_params *p = (struct fn_qaws_params *) params;
  203.   gsl_function *f = p->function;
  204.   gsl_integration_qaws_table *t = p->table;
  205.  
  206.   double factor = 1.0;
  207.   
  208.   if (t->beta != 0.0)
  209.     factor *= pow(p->b - x, t->beta);
  210.  
  211.   if (t->nu == 1)
  212.     factor *= log(p->b - x);
  213.  
  214.   return factor * GSL_FN_EVAL (f, x);
  215. }
  216.  
  217.  
  218. static void
  219. compute_result (const double * r, const double * cheb12, const double * cheb24,
  220.         double * result12, double * result24)
  221. {
  222.   size_t i;
  223.   double res12 = 0;
  224.   double res24 = 0;
  225.   
  226.   for (i = 0; i < 13; i++)
  227.     {
  228.       res12 += r[i] * cheb12[i];
  229.     }
  230.   
  231.   for (i = 0; i < 25; i++)
  232.     {
  233.       res24 += r[i] * cheb24[i];
  234.     }
  235.   
  236.   *result12 = res12;
  237.   *result24 = res24;
  238. }
  239.